function [P,grid_y]=TRP_AR(rho,mean_y,var_eps,no_y, min_y, max_y)
% Purpose: Creates a no x no matrix of discrete conditional transition probabilities
%          assuming AR(1) process of the form
%          y(t)=mean_y+rho*y(t-1)+eps, 
%          eps ~ N(0,var_eps)

% Inputs:    % rho:       persistence of shocks (can be 1)
%            % mean_y:    conditional mean of the process
%            % var_eps    conditional variance of the process
%            % no_y:      number of elements in vector grid_y
%            % min_y, max_y: width of the grid_y (baseline: 3*std_y)
% Output:    % no x no matrix of transition probabilities
%            % no x 1 vector of discrete grid

if nargin < 5
   bandwidth_y = 3; % number of standard deviations around average   
end

if var_eps<0, error('Negative variance not possible!'), end

if rho<1
    std_y=sqrt(var_eps/(1-rho^2));
elseif rho==1
    std_y=sqrt(var_eps);
end

%Grid for the relative prices - defined as (p_i/(A_i*P))
if nargin < 5
   bandwidth_y = 3; % number of standard deviations around average
   min_y    =   mean_y-(bandwidth_y*std_y);
   max_y    =   mean_y+bandwidth_y*std_y;
end
inc_y    =   (max_y-min_y)/(no_y-1);
grid_y   =   (min_y:inc_y:max_y)';

%Calculating conditional weights
y_matr     =   repmat(grid_y,[1 no_y-1]);
mpoints_y  =   repmat((grid_y(1:no_y-1)'+grid_y(2:no_y)')/2,[no_y 1]);

%Cumulative distribution function
cdf_mpoints = normcdf(mpoints_y-rho*y_matr,mean_y*ones(no_y,no_y-1),sqrt(var_eps)*ones(no_y,no_y-1));

%Transition matrix (transposed, to be consistent with the tauchen code)
P = ([cdf_mpoints(:,:) ones(no_y,1)]-[zeros(no_y,1) cdf_mpoints(:,:)])';